In [13]:
import pandas as pd
import os
import yaml
from sentinelsat import SentinelAPI, geojson_to_wkt
import osmnx as ox
In [14]:
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import os
from osgeo import gdal
import rasterio
from rasterio.plot import show_hist
from rasterio.plot import show
from rasterio.mask import mask
from shapely.geometry import box
from fiona.crs import from_epsg
In [15]:
from shapely.geometry import box
In [16]:
with open("auth.yaml", 'r') as ymlfile:
auth = yaml.safe_load(ymlfile)
In [17]:
api = SentinelAPI(auth['copernicus_open_access_hub']['username'],
auth['copernicus_open_access_hub']['password'],
'https://scihub.copernicus.eu/dhus')
In [18]:
# this is the place where I was born and grew up
# I chose to use the image from this place because it would be
# easier to label the clusters since I am closely familiar with the area
place_name='Gelendzhik, Russia'
In [19]:
# will return either polygon or point (in GeodatFrame format) if boundaries of the place aren't available
place = ox.gdf_from_place(place_name)
In [161]:
place
Out[161]:
In [158]:
place.plot()
Out[158]:
In [20]:
place.crs
Out[20]:
In [21]:
query_box=box(minx=place.bbox_west.values[0],
miny=place.bbox_south.values[0],
maxx=place.bbox_east.values[0],
maxy=place.bbox_north.values[0])
In [22]:
# 3-month period to search for scenes with less than 6% cloud coverage where bounding box of the area
# of interest is within the image (that is image contains the bounding box)
products = api.query(query_box.wkt,
area_relation='Contains',
date = ('20190601', '20190926'),
platformname = 'Sentinel-2',
processinglevel = 'Level-2A',
cloudcoverpercentage = (0,5))
In [23]:
products_gdf = api.to_geodataframe(products)
In [24]:
if not products_gdf.empty:
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
else:
print ('No records scences found for the filtered parameters')
In [162]:
products_gdf_sorted.head()
Out[162]:
In [26]:
ax=products_gdf_sorted.plot(column='uuid', cmap=None, alpha=0.5)
place.plot(ax=ax, edgecolor="white", facecolor='None')
Out[26]:
In [27]:
# this is the scene ID needed to to call download fucntion
uuid=products_gdf_sorted.index.to_list()[0]
In [28]:
# this is the name that will be given by default to a folder in which data will be donloaded
title=products_gdf_sorted.loc[uuid].title
In [71]:
api.download(uuid, directory_path='./imagery')
Out[71]:
In [72]:
import zipfile
with zipfile.ZipFile(f"imagery/{title}.zip","r") as zip_ref:
zip_ref.extractall("imagery")
In [144]:
# clean-up
if os.path.exists(f"imagery/{title}.zip"):
os.remove(f"imagery/{title}.zip")
In [29]:
# there is lots of folders that get downloaded, actual images is in the GRANULE/IMG_DATA
sentinel_folder=os.listdir(os.path.join('imagery', f'{title}.SAFE', 'GRANULE'))[0]
In [30]:
sentinel_folder
Out[30]:
In [31]:
os.listdir(os.path.join('imagery', f'{title}.SAFE', 'GRANULE', sentinel_folder))
Out[31]:
In [32]:
# one folder for each of the 13 spectral bands--they are in different resolutions
os.listdir(os.path.join('imagery', f'{title}.SAFE', 'GRANULE', sentinel_folder, 'IMG_DATA'))
Out[32]:
In [33]:
folder=os.path.join('imagery', f'{title}.SAFE', 'GRANULE', sentinel_folder, 'IMG_DATA')
In [34]:
os.listdir(folder)
Out[34]:
In [35]:
for sub_folder in os.listdir(folder):
print (f'{sub_folder} :', )
for file in os.listdir(os.path.join(folder, sub_folder)):
print (file)
Band designations are found here:
In [36]:
resolution='R20m'
In [37]:
os.listdir(os.path.join(folder, resolution))
Out[37]:
In [38]:
red=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B04' in band][0]
blue=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B02' in band][0]
green=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B03' in band][0]
nir=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B8' in band][0]
band5=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B05' in band][0]
band6=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B06' in band][0]
band7=[band for band in os.listdir(os.path.join(folder, resolution)) if 'B07' in band][0]
In [43]:
# read RGB composite to take a look at what we downloaded
src=rasterio.open(os.path.join(folder, resolution, 'T37TDK_20190916T082559_TCI_20m.jp2'))
In [47]:
data=src.read(1)
In [48]:
src.meta
Out[48]:
In [49]:
# make GeodataFrame from the area of interest box
geo_box = gpd.GeoDataFrame({'geometry': query_box}, index=[0], crs=place.crs)
In [55]:
# Vizualize the raster band and the area of interest
fig, ax = plt.subplots(figsize=(12, 12))
plt.imshow(data, extent=[src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top])
geo_box.to_crs(src.crs).plot(ax=plt.gca(), edgecolor="red", facecolor="none");
In [57]:
def crop_band(clip_gdf, input_folder, band, output_folder, out_name):
'''Clip raster band to the area of interest and write out clipped raster (in tif format) to the specified folder
Params:
clip_gdf (GeodatFrame): Area of interest
input_folder (str): Folder with raster band
band (str): Name of the band
output_folder (str): Folder where clipped features will be stored
outname (str): Name for the clipped band
'''
import json
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# this is how the mask function needs the clip feature to be according to the documentation
crop_feature = [json.loads(clip_gdf.to_json())['features'][0]['geometry']]
src=rasterio.open(os.path.join(input_folder, band))
# clip using mask
out_image, out_transform=mask(src, crop_feature, crop=True)
# get the metadata of the original raster to use for the output raster
out_meta=src.meta.copy()
# update the width an height since they are now (after clipping) they are
# different from the original raster's dimensions
out_meta.update({'height': out_image.shape[1],
'width':out_image.shape[2],
'transform': out_transform})
# open the new tif file in the writing mode and write out the clipped raster
with rasterio.open(f'{output_folder}/{out_name}.tif', 'w', **out_meta) as dest:
dest.write(out_image)
In [58]:
band_names={red:'red',
green: 'green',
blue: 'blue',
nir: 'nir',
band5: 'band5',
band6: 'band6',
band7: 'band7',
}
In [60]:
band_names
Out[60]:
In [61]:
# reproject area of interest to same projection as the imagery
geo_box=geo_box.to_crs(src.crs)
In [62]:
os.path.join(folder, resolution)
Out[62]:
In [63]:
for band, name in band_names.items():
crop_band(clip_gdf=geo_box,
input_folder=os.path.join(folder, resolution),
band=band,
output_folder='clipped',
out_name=name)
Some feature generating to add additional dimension¶
Normalized Difference Vegetation Index (NDVI) quantifies vegetation by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs). NDVI always ranges from -1 to +1. But there isn’t a distinct boundary for each type of land cover.
For example, when you have negative values, it’s highly likely that it’s water. On the other hand, if you have a NDVI value close to +1, there’s a high possibility that it’s dense green leaves.
But when NDVI is close to zero, there isn’t green leaves and it could even be an urbanized area.
https://gisgeography.com/ndvi-normalized-difference-vegetation-index/
$NDVI=\frac{NIR-RED}{NIR+RED}$
In [64]:
# write the function to calculate NDVI
def calculate_ndvi(red, nir):
''' Given Red and NIR bands as input,
calcuates NDVI'''
# beacase division will unlikely result in integers
# change the datatype of our arrays from integer to float
red=red.astype(float)
nir=nir.astype(float)
# since we have 0 values
# we need to mask them out, so that they don't participate in the calculations
masked_red=np.ma.masked_where(red==0, red)
masked_nir=np.ma.masked_where(nir==0, nir)
# NDVI formula
ndvi=(masked_nir-masked_red)/(masked_nir+masked_red)
ndvi=ndvi.filled(fill_value=-1)
return ndvi
In [65]:
tifs=[tif for tif in os.listdir('./clipped') if tif.endswith('.tif')]
In [66]:
red=rasterio.open(os.path.join('./clipped', 'red.tif')).read(1)
blue=rasterio.open(os.path.join('./clipped', 'blue.tif')).read(1)
green=rasterio.open(os.path.join('./clipped', 'green.tif')).read(1)
nir=rasterio.open(os.path.join('./clipped', 'nir.tif')).read(1)
band5=rasterio.open(os.path.join('./clipped', 'band5.tif')).read(1)
band6=rasterio.open(os.path.join('./clipped', 'band6.tif')).read(1)
band7=rasterio.open(os.path.join('./clipped', 'band7.tif')).read(1)
In [67]:
red.shape
Out[67]:
In [68]:
ndvi=calculate_ndvi(red, nir)
In [79]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
In [81]:
fig, ax = plt.subplots(figsize=(12, 12))
im=plt.imshow(ndvi, cmap=plt.get_cmap('RdYlGn'))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(im, cax=cax);
In [100]:
from sklearn import cluster
In [101]:
all_bands=np.dstack([red.ravel(),
green.ravel(),
blue.ravel(),
nir.ravel(),
ndvi.ravel(),
band5.ravel(),
band6.ravel(),
band7.ravel()])
In [102]:
all_bands.shape[1:]
Out[102]:
In [103]:
X = all_bands.reshape(all_bands.shape[1:])
In [104]:
X.shape
Out[104]:
In [105]:
k_means = cluster.KMeans(n_clusters=6, random_state=123)
k_means.fit(X)
Out[105]:
In [88]:
X_cluster = k_means.labels_
In [89]:
set(k_means.labels_)
Out[89]:
In [90]:
X_cluster.shape
Out[90]:
In [91]:
X_cluster=X_cluster.reshape(red.shape)
In [94]:
plt.figure(figsize=(12, 12))
ax = plt.gca()
im=plt.imshow(X_cluster, cmap=plt.cm.get_cmap('cubehelix', 6))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ticks=range(6), label='cluster', cax=cax)
plt.show()
In [96]:
from PIL import Image
In [97]:
# I used GoogleEarth to label each cluster
# I zoomed into each area for each cluster to give a meaningful cluster label
image=Image.open(os.path.join('./clipped', 'GoogleEarth.png'))
In [99]:
fig, ax = plt.subplots(figsize=(11, 11))
ax.imshow(image, interpolation='nearest')
plt.tight_layout()
In [106]:
labels={
3: 'trees/dense vegetation',
1: 'water',
0: 'residential/roads',
4: 'grassland',
2: 'commercial/industrial',
5: 'sparse vegetaion' ,
}
In [107]:
bands_df=pd.DataFrame(X, columns=['red', 'green', 'blue', 'nir', 'ndvi', 'band5', 'band6', 'band7'])
In [108]:
bands_df['cluster'] = k_means.labels_
In [109]:
bands_df['labels']=bands_df['cluster'].map(labels)
In [110]:
bands_df.sample(20)
Out[110]:
In [111]:
import seaborn as sns
In [112]:
sns.pairplot(bands_df.drop('cluster', 1), hue="labels")
plt.show()
In [113]:
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
In [114]:
X_bands=bands_df.drop(['labels', 'cluster'], 1)
In [115]:
y=bands_df.cluster
In [116]:
ss=StandardScaler()
In [117]:
Xs = ss.fit_transform(X_bands)
In [119]:
X_train, X_test, y_train, y_test=train_test_split(Xs, y)
In [120]:
grid_params={'n_neighbors': list(range(3, 50, 4)),
'weights': ['uniform', 'distance'],
'metric': ['euclidean', 'manhattan']
}
In [121]:
gs=GridSearchCV(KNeighborsClassifier(),
grid_params,
verbose=1,
cv=3,
n_jobs=-1
)
In [122]:
gs_results=gs.fit(X_train, y_train)
In [137]:
gs_results.best_score_
Out[137]:
In [138]:
gs_results.best_params_
Out[138]:
In [139]:
# now train the model using the best parameters returned by Gridsearch
optimal_knn=KNeighborsClassifier(n_neighbors=gs_results.best_params_['n_neighbors'],
weights=gs_results.best_params_['weights'],
metric=gs_results.best_params_['metric'])
In [140]:
optimal_knn.fit(X_train, y_train)
Out[140]:
In [141]:
y_pred=optimal_knn.predict(X_test)
In [142]:
from sklearn import metrics
In [143]:
# baseline accurace
y.value_counts() / y.count()
Out[143]:
In [146]:
metrics.accuracy_score(y_test, y_pred)
Out[146]:
In [151]:
import joblib
import pickle
In [131]:
joblib.dump(optimal_knn, 'KNeighborsClassifier_LandCover.pkl')
Out[131]:
In [132]:
pwd
Out[132]:
In [133]:
ls
In [157]:
# knn_LandCover = joblib.load('KNeighborsClassifier_LandCover.pkl')
In [ ]: